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ABSTRACT 


A COMPUTER MODEL OF A PHASE -LOCK LOOP 

BY 

RALPH PAUL SHELTON, B.S. E.E. 


Master of Science in Electrical Engineering 
New Mexico State University 
Las Cruces, New Mexico, 1970 
Dr. Frank F. Carden, Chairman 


This thesis presents a computer model of a PLL (phase-lock 
loop), preceded by a bandpass filter, which is valid when the band- 
width of the bandpass filter is of the same order of magnitude as 
the natural frequency of the PLL. This model is shown to agree 
with previously known results for the situation in which PLL 
natural frequency is much smaller than the bandpass filter band- 
width. New results for the PLL natural frequency equal to the 
bandpass filter bandwidth are presented. 

The model presented herein is for a second order PLL operating 
with carrier plus noise as the input. However, it is shown that 
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extensions to higher order loops , and to the case of a modulated 
carrier are straightforward. 

The new results presented give the cycle skipping rate of the 
PLL as a function of the input carrier to noise ratio when the PLL 
natural frequency is equal to the bandpass filter bandwidth. Pre- 
liminary results showing the variation of the output noise power 
and cycle skipping rates of the PLL as a function of the loop damp- 
ing ratio for the PLL natural frequency equal to the bandpass filter 
bandwidth are also presented. 
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CHAPTER I 


INTRODUCTION 

The phase-lock loop is a device with a wide variety of applica- 
tions in electronic systems. F. M. Gardner lists five representative 
applications in his book on the subject [l]. These are? a coherent 
transponder, frequency modulation demodulator, stabilizing oscilla- 
tors, frequency multipliers and dividers, and pulse code modulation 
demodulators. In general, a phase-lock loop is useful whenever one 
needs a noiseless version of a noisy signal which is phase coherent 
with the original signal. 

There is a large body of information concerning many specific 
uses of the PLL (phase-lock loop). However, no general solution, 
that is, one describing the function of the PLL for all situations, 
has been found for a loop of order higher than one. This is a 
direct result of the difficulties involved in solving the nonlinear 
differential equations describing these higher order loops. In 
spite of this difficulty, it is known that a second order PLL is 
superior to a first order PLL for many applications. Enough experi- 
mental data for second order loops has been obtained to allow such 
loops to function better than other known devices for some applica- 
tions . One of the most glamorous applications of second order PLL 
has been in the communications systems for the Apollo program. 

In the applications previously mentioned, the PLL is generally 
preceded by a bandpass filter. This bandpass filter is frequently 
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the IF (intermediate frequency amplifier) of a communications 
receiver. Thus, the. input to the PLL is signal plus band- 
limited no i se. This is the situation modeled in this thesis. 
Whenever the bandwidth of the IF is wide compared to the bandwidth 
of the PLL, one can ignore the bandpass characteristics of the IF 
and analyze the PLL as if the input were white noise plus the 
desired signal. This simplifies the analytical problem consider- 
ably. Most previous work with the second order PLL has been done 
for this case. An important computer model for the white noise case 
has been developed by Sanneman and Rowbotham [2]. This model was 
used to obtain the mean time to unlock for a carrier tracking PLL. 

It will be shown that the results of the computer model developed 
in this thesis agree with this mean time to unlock in the appro- 
priate region. The model presented in this thesis is an improvement 
over previous work because it is valid when the input to the PLL 
cannot be considered to be white noise plus signal. 

The model development begins with the derivation of a general 
baseband model for a PLL. This development, given in Chapter II, 
parallels that of Viterbi [3], Chapter III specializes the model of 
Chapter II to the case with no noise. Computer simulation results 
which verify this specialization of the model are presented. 

Chapter IV is devoted to the development of a digital computer 
representation of the noise process used in the general development 
of Chapter II. Computer results showing the validity of this model 
are presented. Chapter V combines the results of Chapters III and IV 
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to give a computer implementation of the general model of Chapter IT 
for the case ■when no modulation is present. The results of this 
model are compared to previously published work and some recently 
obtained experimental data to show the validity of the model as a 
whole. Chapter V also presents new results obtained with this model. 
Chapter VI gives conclusions and recommendations for future studies 
using this model and extensions thereof. 



CHAPTER II 


DEVELOPMENT! OF THE GENERAL BASEBAND MODEL 

2.1 Introduction 

This chapter develops a baseband model for a phase -lock loop 
preceded by an IF. The development begins with the representation 
of a PLL in terms of ideal physical components. This representation, 
which is standard in the literature, is shown in figure 1 [l,2,3,7, 
15]. The development of the baseband model parallels that of 
Viterbi [3].. The second section of this chapter specializes the 
model slightly. This specialization will simplify the explanation 
of the computer implementation given in Chapters III and V. 

2 .2 General Development 

The development of the general model will begin with the 
mathematical models for the physical devices shown in figure 1. 

The VCO (voltage controlled oscillator) is an oscillator for 
which the output frequency is directly proportional to the input 
voltage. In terms of instantaneous quantities: 

0 Q (t) = K^Ct) (2.2.1) 

where v is the input voltage and 0 Q (t) is the phase of the 
output sinusoid. 

The multiplier is modeled as an approximately ideal device, 
the departure from ideal being that the output is proportional to 
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the product of the inputs. That is, for inputs 1^ and I 2 , and 
output Y , we have: 

Y = K 2 X 1 I 2 „ (2.2.2) 

K 2 is the multiplier gain. Since the multiplier is used as a 
phase detector in this system, K 2 is often referred to as the 
phase detector gain. 

The linear filter is characterized by its impulse response, 
h(t) . In the work following, it will be specialized to a second 
order filter. For the present, it will be most convenient to retain 
the general representation. 

As mentioned in the Introduction, it is desired to model the 
PLL when it is preceded by an IF. In this case, the input may be 
represented as the sum of a signal plus noise. 

x(t) = s(t) + n(t) . (2.2.3) 

The signal considered is a general angle modulated sinusoidal 
carrier . 

s(t) = nT 2A sin [m Q t + 0 i (t)] . (2.2.4) 

The noise is a narrow -band process because of the 
characteristics of the IF. S. 0. Rice [4] has shown that such 
a noise process may be represented as: 

n(t ) - n g (t)sino3 0 t + n c (t)coso) 0 t (2.2.5) 
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where n (t) and n (t) are independent stochastic processes with 
s c 

baseband power spectra identical to that of the original process 
[9], 

The output of the VCO is taken to be; 

z(t) = , '/2B cos [o}Qt+0Q(t)] . (2.2.6) 

Now R(t) can be calculated as follows: 

R(t) = x(t)z(t)K2 • (2.2.7) 

By substitution (2.2.7) becomes: 

R(t ) - K 2 [ \I~2A sin[iD o t+0 i (t)3+n g (t)sina5 o t+ n c (t)coso) 0 t ] 

[ \[ 2 B COSO^t + 0Q (c)] . (2.2.8) 

This result is now simplified using the following 
trigonometric identities: 

sin(x)cos(y) =%sin(x+y) +%sin(x+y) (2.2.9) 

cos(x)cos(y) =%cos(x + y) +%cos(?c+y) . (2.2.10) 


After applying the identities^ (2.2.8) becomes: 
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R(t) = 2ABK 9 [%sin{a} 0 t. + ^(t) + 0J Q t + 0 o (t)} 

+ %s infant + 0 ± (t ) - m Q t - 0 Q (t)}] 

+ NT2BK 2 n s (t)(^sin[coQt +co Q t + 0 Q (t>] + ^sin[a) Q t - a) Q t - ( t ) ] } 

+ ^BK 2 n c (t)(%cos [a)gt + o)gt + 0g(t)] + %cos [a^t - a) Q t - 0 Q (t) ] } . (2 .2 .11) 

Next y(t) is to be calculated. The calculation can be 
simplified by observing at the outset that the linear filter 
is a lowpass filter with a corner frequency much less than oDq . 

This allows one to ignore the terms centered about 2 c0q and 
writes 

R(t) = K 2 ABsin[0 i (t) - 0 Q (t)] + V%K 2 Bn g (t )sin[ “0 Q (t ) ] 

+ */%K 2 Bn c (t)cos [ -0g(t) ] . (2.2,12) 

With R(t) as the input to a linear filter with 
impulse response, h(t), the output may be written as; 


y(t) = / R(u)h(t - u)du 


Now substituting y(t) for v in (2.2.1), 


0 O = K x y(t) 

y(t> = e 0 /\ . 


(2.2.13) 

(2.2.14) 

(2.2.15) 
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Combining (2.2.12), (2.2.13), and (2.2.15) yields i 
e Q - K^K^AB sin[S_^(u) - 0^(u) ]h(t - u)du 

+ n g (u)sin[ -0 q(u) ]h(t - u)du 

+ ^K.tK 0 B / n (u)cos [ -0 n (u) ]h(t - u)du . (2.2.16) 

1 z Jo c u 

This equation is represented by the system shown, in figure 2. 

2 . 3 Some Specializations for Simplification 

There are some specializations of the model, shown in figure 3, 
which simplify the implementation and retain much of the generality. 
First, the product BK^K^ is a combination of loop parameters 
which may be considered as one parameter. Thus, the loop gain K 
is defined as; 

K = BK x K 2 . (2.3.1) 

Second, note that the sine is an odd function and the 
cosine is an even function. Thus; 

sin(-0 Q ) - -sin0Q (2.3.2) 

cos(-e Q ) = cose 0 . (2.3.3) 

Further if one considers only zero mean gauss ian processes 
for n(t), it is immaterial whether one writes n(t) or -n(t) . 
Thus, both inverters in this portion of the loop may be omitted. 
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With the two previous simplifications one may model the system 
as shown in figure 3. The equation describing the system then 
becomes : 


0 Q = AK^ sin[0^(u) - £q(u) ]h(t - u)du 


n g (u)sin0Q(u)h(t - u)du 


+ V *jf n c (u)cos0 o (u)h(t - u)du . (2.3.4) 




The PLL in Terms of Physical Components 
figure 1 



Block Diagram of the General Baseband Model 
figure 2 



Simplified Block Diagram of the Baseband Model 

figure 3 













CHAPTER III 


DETERMINISTIC MODEL VERIFICATION 

3 . 1 Introduction 

In this chapter the general model of Chapter II will be 
specialized to a second order phase- lock loop. Further, only the 
deterministic response will be considered. The implementation of 
the deterministic model will be discussed in detail. Finally, the 
results obtained with this implementation will be shown to agree 
with previously published results for the deterministic response of 
a second order PLL , 

3 . 2 Deterministic Model 

A model for a phase- lock loop with deterministic inputs may be 

derived, beginning with (2.3.4). If n (t) and n (t) are set to 

s c 

zero in (2.3.4) one obtains; 

% AK Jq sin ^ i ( u ) " 0 o ( u )]h(t - u)du . (3.2.1) 

This equation is represented in block diagram form in figure 4. 

Figure 4 may also be viewed as being derived directly from figure 3 

with n (t) and n (t) set to zero, 
s c 

This is a standard representation for the deterministic base- 
band model, and similar developments may be found in books by 
Viterbi [3] and Gardner [l]. 
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All results presented in this thesis are based on the following 
specializations. First, the carrier amplitude A is assumed fixed 
at unity. The filter, represented in general by h(t) in (3.2.1), 
was specialized to a filter of the form* 

H,. (s) = 1 + - (3.2.2) 

t s 

in terms of the s- domain transfer function. 

Now rewriting (3.2.1) in s -domain notation, substituting 
(3,2,2) for H(s), and setting A = 1, one obtains: 

se Q (s) * K[sin{e i (s) - @q(s)} ] [l + . (3.2.3) 

Thus, (3.2.3) is represented in block diagram form in figure 5. 

Figure 5 represents the way in which the deterministic PLL 

model was implemented. A few more steps are necessary to represent 

this model in terms of the parameters most commonly used in the 

literature. The first step in this procedure will be to assume 

0. » 0 A and hence sin(0. - 0 A ) ~ 9 . - 0 A . Thus, (3.2.3) now 

iO i 0 i 0 7 

becomes: 

s0 Q (s) = K[0.(s) - 0 Q (s) ] [1 + ~] , (3.2.4) 


This may be rewritten as 


e i (s ) 


sK + aK 
s 2 + Ks + aK 


(3.2.5) 
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Transfer functions with second order denominators commonly 
occur in elementary closed loop systems theory. The response of 
such systems is characterized by two parameters . These are normally 
referred to as zeta (£), the damping ratio, and to , the system 
natural frequency. The denominator of the transfer function is 
normally written as: 

s 2 + 2Ca> s + to 2 . (3.2.6) 

5 n n 

Thus, for the second order PLL: 


K = 2£tjo 


(3.2.7) 



(3.2.8) 


In terms of these parameters, (3.2.4) can be rewritten 
as : 


s6 0 (s) 


[e i (s) - e Q (s)] 




(3.2.9) 


and figure 5 can be redrawn as shown in figure 6. 

The form of the model depicted in figure 6 was used for deter- 
ministic studies of the second order PLL. 
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3 , 3 Computer Implementation 

It will be most convenient to describe the digital computer 
implementation of the model with reference to a block diagram., 

Figure 6 represents the system simulated. Since the simulation is 
in the time -domain, figure 6 is redrawn, using time -domain notation, 
as figure 7. 

The digital computer simulation was derived directly from this 
block diagram. The operations indicated on the diagram were per- 
formed over discrete blocks of time. A computer program for step 
response and a block diagram similar to figure 7 annotated with the 
computer program variables are included in Appendix A. 

Three important points with respect to the simulation will now 
be considered. First, the basic simulation rate will be considered; 
second, the integration method; and third, the determination of the 
input . 

The basic simulation rate was empirically determined by starting 
at a very high sample rate, with a step input to the system, and then 
lowering the sample rate until some degradation was observed. In 
this context degradation means a change in the output was observed. 
This is based on the assumption that as one increases the sample 
rate a better representation of the system is obtained. The sample 
rate was then set. at a slightly higher rate. This rate may also be 
viewed as the integration interval for it is the basic time segment 
for the numerical, integration technique. For all results presented 
in this thesis, the sample rate was 20 samples per second. 
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s -domain Block Diagram of Deterministic Model 

figure 6 



Time-domain Block Diagram of Deterministic Model 

figure 7 
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The integration method used was suggested by F, E. Nixon [5] . 

It was chosen simply because it is easy to implement. This method, 
called the weighted average method, averages the derivative of the 
integral over two successive time intervals and uses this value as 
the rate of change of the integral. The use of an average deriva- 
tive is thought to add stability to the numerical technique. How- 
ever, no detailed analysis of the technique was made. The 
integration method is not. claimed to be optimal but merely sufficient 
for the needs of the work presented here. However, some preliminary 
comparisons indicate this method is much faster on the digital com- 
puter than more commonly used differential equation solution methods 
such as the Runge-Kutta method. A study to determine which of the 
many available numerical integration techniques is optimum in some 
particular sense could be a topic worthy of a thesis in itself. 

The input to the model just developed is the phase of the 
input carrier. It should be noted that references to the type of 
modulation refer to frequency modulation. Thus, a step input is a 
step in frequency which is a ramp in phase, and a cosine modulation 
appears at 0^(t) as a sine. This is an important point to remem- 
ber when interpreting the programs included. 

3 .4 Verification of Deterministic Response 

One of the primary areas of interest in this investigation is 
the cycle skipping behavior of the phase- lock loop. Thus, the 
investigations for deterministic inputs were designed to show that 
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the phase- locked loop, as modeled in this presentation, loses lock 
for the same deterministic inputs as in previous studies. 

The results for a frequency modulation step are shown in 
figure 8, The heavy lines are from reference [7], The points 
shown with circles were determined with the model described earlier 
in this chapter . 

Figure 9 shows regions of stability and instability as presented 
in reference [6] for frequency modulation by a single sine wave. 
Points determined using the model of this chapter are also shown in 
figure 9. Circles represent stable points and triangles represent 
unstable points. 

The agreement shown by the figures is very good. It should be 
noted that reference [6] mentioned that the line shown in figure 9 
does not represent a sharply-defined boundary. Further, the center 
point in each group of three points had rather large phase errors 
(greater than -) . However, the output was still a reasonable 
approximation of a sine wave. 
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Damping Ratio (£) 


Verification of Deterministic Model for Step Input 

figure 8 



Verification of Deterministic Model for Sine Input 

figure 9 


CHAPTER IV 


THE NOISE MODEL 


4 «1 Introduction 

In this chapter a model for the narrow -band noise processes, 

n (t) and n (t) of Chapter II, will be developed . The model is 
s c 

specifically designed for easy implementation on a digital computer. 
The model begins with a wide-band noise source which is described in 
section 4.2. The wide-band noise from this source is passed through 
a filter with a passband narrow compared to the spectrum of the 
noise source. Thus the output of the filter is the same as if the 
input were white noise. The filter is described in section 4.3. 

The results of tests on the noise model as a whole are presented in 
section 4 .4 . 

4 „ 2 Th e Wide-band Noise Source 

The fundamental noise source was a gauss ian random number 
generator. The output numbers from this source were independent 
with a gaussian probability distribution about a zero mean. This 
random number generator is described in Appendix B. Each number 
from the random number generator was taken to represent the ampli- 
tude of a time function for a period of At seconds. This time 
function is illustrated in figure 10. 

Now the power spectrum of the process will be found. This is 
the Fourier transform of the autocorrelation function. For the 
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following calculations the process will be assumed ergodic . The 
ergodicity of the process will be used to calculate the autocorrela- 
tion function. This development follows that of Lee [8]. 

For an ergodic process , x(t), the autocorrelation function may 
be written as : 

R(t) = lim f x(t)x(t. + t)dt. . (4.2.1) 

T-m» 


This integral will be evaluated over two regions and the results 
combined to give R(t) for all t . 

Firsts consider the region where ] t | > At . For this region, 
x(t) and x(t + t) are independent. Since the process is 
generated with zero mean, the average of the product of x(t) and 
x(t + x) will be zero for this region. Hence: 

R(t) = 0 | t |>A t . (4.2.2) 

To evaluate R(x) for |x|<At it will be most advantageous 
to first consider a single region of width At. The functions 
x(t) and x(t + t) are illustrated for a typical region of this 
type in figure 11. The number is the output of the random 

number generator in the interval At . From the illustration the 
integral of the desired product in the interval At can be viewed 
as an area and written as: 


x(t)x(t+ x)dt = N.N. _t + N. (At - t) 
x i-l 1 v 


1 


h+ At 


(4.2.3) 
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Over this interval the normalizing factor multiplying the integral 
is simply “ . Hence, for the interval, the integral -is; 

^ N i N i-l T ]/ At + ^ N t 2 ^ At T )^/ At • (4.2.4) 

The entire integral corresponds to averaging over all such 
intervals. The average of N.N. - is zero from the way in which 

1 It 1 

2 

the N^ are generated. Further, the average of FL is simply 
2 

the variance, o , of the random number generator. Hence, applying 
the same arguments to negative shifts : 

R(t) - GT 2 (At - | t | ) / At | t |<A t . (4.2,5) 

This result is shown pictorially in figure 12. 

Now for the power spectrum, S^Cct), the Fourier transform 
of R(t) is taken 


S (<x>) = / R(t) e ^^dx 

n J- CO 


(4.2.6) 


for this R(x): 


S (o>) = 
n 


r° g 2 (At + t) e -ja)x dT+ 

JAt At 


'At 


e^dT .(4.2.7) 


This transform is a well-known result. (See for example 
reference [9] page 340.) For S n (a}) one obtains: 

s (o>) = 

" oAt 
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This result is illustrated in figure 13. 

The computer model was implemented with At = .05 and 

2 

a - 1.0 . This gives the power spectrum a peak amplitude of .05 
and the first zero crossing at 125.5 radians. 

4 .3 The Lowpass Fil ter 

The lowpass filter chosen was a fifth order Ghebyshev filter 
with a 0,5 dB passband ripple. The 0.5 dB cutoff frequency was set 
at 1.0 radian. Thus, the power spectrum of the noise process 
previously described is essentially flat (the variation is less 
than one percent) with an amplitude of .05 across the passband of 
the filter . 

The implementation of the filter was similar to an implementa- 
tion which could be used on an analog computer. That is, the 
differential equation describing the filter was solved using five 
integrators in series. The transfer function realized, in terms of 
s-domain notation, was: 

H C (S)=^ r J= 2 • (4.3.1) 

s + b. s + b 0 s + b^,s + b_ s +b_ 

4 3 2 1 0 

A block diagram of the method used is shown in figure 14. The 
coefficients, , were taken from Weinberg [10], A listing 
of the computer program generating the filtered noise is included 
as Appendix C . 

The frequency domain characteristics of the filter, as 
implemented on the digital computer, were determined by passing 



Wide -band Noise Power Spectrum 
figure 13 


C> |to 
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sinusoids of known amplitude through the filter. The results of 
this test were found to be very close to that of an ideal Chebyshev 
filter. These results are plotted in figure 15. The plot is for a 
normalized filter. 

There were many reasons for choosing this particular filter. 
The integration technique used was the same as that described in 
section 3.3. This was the highest order filter which could be 
easily implemented on the digital computer with this simple inte- 
gration technique. It is sufficiently close to an ideal square 
filter to allow many analytical calculations to be done with this 
approximation. It is also very close to the characteristics of 
many IFs found in practical systems. 

4 .4 Empirical Verification of the Noise Model 

This section describes tests made on the entire noise model 
which indicate that the previous analytical work; including 
assumptions, is valid. 

The first test run was an estimation of the mean and variance 
of the output noise. This was accomplished by calculating the mean 
and variance of samples of the noise generated. For this test the 
output of the filter was multiplied by b^ to give a normalized 
result . 

The input process was generated with zero mean. Thus, the 

output of the linear filter should have zero mean. This was found 

to be consistently true, to within the accuracy of the digital 

5 

computer, for samples of 2 x 10 points. 
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The variance of the output process can be calculated by assuming 
the input spectrum to be flat. 

r°° 

^ J-ca s 0 (a))d ^ • (4,4.1) 

Sq(oO is the power spectrum of the output process . For a 
linear system; 

Sq(o>) - f H( joo> | 2 S i (o}) . (4.4,2) 

Thus, for this system; 

p 00 

V - %* / N Q /2 |H( jcu) 1 2 d ui . (4.4.3) 

With the assumption of an ideal square filter of width one 
radian per second, the output variance becomes: 

o Q 2 ~ n q / 2]-2 . ( 4 . 4 , 4 ) 

From the results of section 4.2; 

a Q 2 = .01590 . (4.4.5) 

The integral: 





|H(co) | 2 dcu 


(4.4.6) 


has been evaluated for a fifth order Chebyshev filter [ll] . 
Using the tabulated value; 
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ct q 2 = .01625 


(4.4.7) 


averaged over 4 x 10 sample points it was found that: 


a Q 2 = .01610 


( 4 .4 . 8 ) 


This is very good agreement. 

Since the determination of the rate of impulses produced by a 
phase-lock loop tracking a carrier plus noise is an important 
objective of this work, the rate of impulses in the IF noise 
model is important. S. 0. Rice [4] has developed a formula 
for this rate based on the radius of gyration of the noise: 

- r(l- erf n / p) impulses per second (4.4,9) 



where p is the numerical carrier to noise ratio and r is 
the radius of gyration of the noise. The erf(x) is defined 
by: 


erf (x) 



(4.4.10) 


In order to better explain this calculation, a short 
digression on the definition of an IF impulse is in order. 

A phasor representation of the carrier plus noise is shown in 
figure 16. According to Rice's definition, an IF impulse occurs 
whenever the resultant phasor crosses the zero axis in the left 
half plane. The carrier amplitude in this model is set to V2 . 
Thus, Rice's impulses are counted with this carrier amplitude. 



Normalized Amplitude (dB) 
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Rice gives the radius of gyration of noise with an ideal square 
power spectrum as: 

r = £/ \fl2 (4. 4 ..11) 

where (3 is the IF bandwidth. For an ideal square filter with 
a cutoff of one radian per second: 

r ± = .0902 . (4.4.12) 

If the spectrum of the noise is assumed to be that produced 
by an ideal fifth order Chebyshev filter, the radius of gyration 
can be calculated: 


r = .1525 . (4.4.13) 

c 

This calculation is performed in Appendix D. 

An empirical determination of the radius of gyration of the 
noise can be made by taking one-half the zero crossing rate of the 
quadrature component of the noise, [4]. This estimate for the 
radius of gyration was made by observing the component for 4 x 10 
seconds and gives an empirical radius of gyration of; 

r * .1066 . (4.4.14) 

e 

This indicates the variation of the filter, as implemented, from 
the ideal Chebyshev filter. 

These rates and the rate observed during simulations are plotted 
versus carrier to noise ratio in figure 17. The IF impulse rate 
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observed during simulation is almost exactly that predicted with the 
empirically determined radius of gyration. Further the rate pre- 
dicted by the empirically determined radius of gyration lies above 
that predicted by an ideal square filter and below that predicted 
for the ideal Chebyshev filter. 



Impulses per Second 
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CHAPTER V 


THE GENERAL COMPUTER MODEL 


5 .1 Introduction 

This chapter gives a detailed description of the digital com- 
puter implementation of the general PLL model developed in Chapter 
II. The method of inserting noise in the model is described. The 
method of determining the occurrence of PLL impulses and IF impulses 
is described . Section 5.4 presents cycle skipping rates (impulse 
rates) for the case in which the IF bandwidth is much greater than 
the PLL bandwidth and compares them to previously published 
results . 

The relationship between loop carrier to noise ratio and output 
noise power is compared to recently obtained experimental results. 
Section 5.5 presents some new results obtained with the model, pre- 
sented in this thesis. These illustrate ways in which the model may 
be used. 

The model was implemented in two phases. First, two time series 
of noise were generated on magnetic tape. Second, these noise pro- 
cesses were read into the PLL modeling program and used for n (t) 
and n c (t) which are described in Chapter II. This approach was 
necessary because large amounts of time were required to generate 
the noise. A program was written which generated the noise it used 
internally. This approach was over 12 times as slow as the simula- 
tion using the noise from tape. 
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The noise processes, n (t) and n (t), are assumed independent 

s c 

in the derivation of the model. This vas guaranteed in the compu- 
ter model by the fact that the processes were generated in blocks 
of 2 x 10^ points (see Appendix B) . The last random number in the 
random number generator was saved at the end of each run and used 
to initialize the random number generator for the next run. The 
random numbers were independent of one another and thus the only 
correlation was that produced by the filter. The 2 x IQ"* points 

4 

correspond to 10 seconds which is many times the decorrelation 
time of the filter. Thus, corresponding points in all of the time 
series are uncorrelated. 

The same time series of noise from tape were used for many 
runs. The variations in carrier to noise ratio were obtained by 
multiplying the noise from tape by appropriate gain factors. The 
calculation of the gain factor is described in detail in section 5.2. 

While the simulation was running, the test described in sec- 
tion 5.3 determined the occurrence of impulses in the IF noise. 
Section 5.4 describes the method used to detect PLL impulses. The 
time of occurrence of each IF impulse was printed as was the time 
of occurrence of each PLL impulse. The duration of each IF impulse 
was also part of the output . Finally, the total number of impulses 
of both types was summarized at the end of each run. The average 
noise power out of the PLL was also calculated for each run. The 
output format was such that additional information may be gleaned 
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from the data collected. This information is not directly appli- 
cable to the specific objective of this thesis . 

A complete listing of the main simulation program is included 
in Appendix E . 

5 . 2 The Gain Factor 

A gain factor (called GAIN in the FORTRAN program) was used to 
adjust the amplitude of the noise time series to give the carrier to 
noise ratio desired in each simulation run. The components of this 
factor will be presented in three parts. First,, there is an adjust- 
ment for the Chebyshev filter gain; second is the main portion of 
the factor which includes the carrier to noise ratio, and finally a 
factor of v 1/2 introduced by the way in which the noise appears in 
the baseband model . 

In the developmental stages of this work it was found that the 
simplest way to view the filter was as a normalized filter followed 
by an ideal amplifier. The noise generated on tape was the output 
of an unnormalized filter with gain, l/b Q . Thus, the first portion 
of the gain factor was b^ which effectively normalized the filter. 

In section 3.2, the model was specialized by setting the rms 
values of the carrier amplitude to one. Thus, the CNR (carrier to 
noise ratio in the IF) is simply the reciprocal of NP (noise power 
out of the fifth order Chebyshev filter) . 


CNR = 1/NP . 


(5.2.1) 
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The wide-band noise process has been shown to be flat across 
the passband of the filter with an amplitude in the passband of At . 
The noise power out of the filter is: 

NP 0 = 2n : ^ oo S 0^ a ^ dm * (5,2.2) 


For a linear system with an input of white noise with spectral 
density, Nq/ 2, the output power spectrum is given by: 

S Q to) = [N 0 /2] | H( ja>) | 2 . (5.2.3) 

Thus, substituting (5.2.3) in (5.2.2) 

i r N o 

NP 0 = T H(jcu)H( ' j£u)d “ • (5.2.4) 


Now a change of variables, substituting s for jco , gives: 


NP Q = N q /2 


,+joo 


2*j 


hy 


H(s)H( -s)ds 


(5.2.5) 


where H(s) is the s-domain transfer function of the filter. 

The integral in brackets (defined I for convenience) has 

m 

been tabulated for a fifth order Chebyshev filter [ 1 1 ] . 

I m = 1.02/* . (5.2.6) 

For convenience the approximation 1.02 ” 1.0 will be used. With 
this approximation, substituting (5.2.6) in (5.2.5) gives: 
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NP Q = N 0 /2* . (5.2.7) 

Next, note that passing the output of the filter through an 

ideal amplifier with gain, G , will multiply the output power 
. 2 

hy G . Using such an amplifier in the model to provide 
the necessary amplitude changes gives: 

NP = G 2 N q /2* . (5.2.8) 

Substituting (5.2.8) in (5.2.1) and solving for G gives: 

G - . (5.2.9) 

This is the gain factor which would be needed on the output 
of a normalized filter to give a particular CNR. 

Finally, from Chapter II, it can be seen that when the 
IF noise processes are n (t) and n (t), the noise processes 

S G 

in the phase-lock loop model are n (t) and \[% n (t) . 

s c 

Thus, the gain factor for the PLL model must be divided by the *7*2 . 

Combining the three factors just described gives the gain 
factor found in the PLL modeling program. 

GAIN = .178923 \T 20jt/ 2CNR . (5.2.10) 

In (5.2.10), .178923 is b^ for the fifth order Chebyshev 
filter [10]. 
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5 .3 IF Impulse Detector 

The definition of IF impulses has been described in detail in 
section 4.4. The most significant point to be brought out here is 
related to the gain factor previously calculated. The noise pro- 
cesses of the IF were divided by the v2 for input to the PLL 
model. It is much faster, in terms of computer time, to determine 
the IF impulses with the scaled version of the noise than to rescale 
the noise for this operation. 

From the phasor diagram of figure 16 one can see that no IF 
impulse can occur while the amplitude of n^(t) is greater than the 
negative of the carrier amplitude. This condition was sensed by 
logic in the program. When this condition existed, each sign change 
of n g (t) was counted as an IF impulse. 

In section 2.3 the carrier amplitude was fixed at J~2 . Thus, 
for an IF impulse to occur, the condition: 

n.(t) < - 41 (5.3.1) 

must exist. This condition implies that, in terms of the 
input noise processes for the PLL: 

n c (t )/ *J1 < -1 . (5.3.2) 

This was the test used in the program logic. The sign changes of 
n g (t) are the same as those of n g (t) . Thus, when the first 
condition was satisfied, the sign of the product of the current and 
previous points of n g (t) was tested. If the sign was negative, 



an IF impulse was said to have occurred. One should note that the 
definition of an IF impulse used here is that given by Rice [4]. 
Other definitions have been used and could yield different results. 

5 .4 PLL Impulse Detector 

PLL impulses were detected by observing a short term average of 
the output phase of the FLL . When comparing the results of this 
thesis with other work in the area, one should note that the output 
phase in this model is the negative of the difference in the phase 
of the VCO and the phase of the carrier. This result is obtained 
because carrier phase is set to zero in this model; that is, no 
modulation is present. 

The differential equation describing the PLL, for no modulation 
present, has been shown to have stable points for an output phase of 
± 2n* where n = 0,1, 2,3 . . . . [3] . In the case with no modulation, 
the PLL produces an output impulse when the PLL moves from the 
vicinity of one stable point to another. This event requires a 
change in the average value of the output phase of 2 Jt . The output 
of the PLL is the derivative of the output phase, and hence a level 
change in the output phase corresponds to an impulse at the PLL 
output . 

In the early stages of the development of the model, the output 
phase was printed at each simulation interval. It was observed that 
whenever the long term average of the output phase changed by 2rt , 
the change was very definite. Furthermore, the initial transient 
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occurring during the move from one stable point to another carried 
the output phase beyond the next stable point* 

In the model presented in this thesis , the output phase was 
averaged over intervals two seconds long* Whenever this average was 
greater than 2^ from the stable point around which the loop had been 
operating, the loop was said to have moved to a new stable point* 

This new value of ± 2nrr was used for future comparisons * The two 
second averaging interval was chosen on the basis of the fact that, 
when the output phase was continuously printed, the output phase 
never went to a new stable point and stayed in that vicinity for this 
length of time and immediately returned to the original stable point. 
Qualitative analysis of the output phase points indicated that the 
changes from one stable point to another were sufficiently definite 
that the choice of the averaging interval was not critical, 

5 . 5 Comparison with Other Results 

The first results presented for comparison were originally pre- 
sented in a paper by Sanneman and Rowbotham [2], Through a computer 
simulation they determined the mean time to unlock for a PLL operat- 
ing with carrier plus white noise as the input. Their result is 
given as a function of the CNR^ (carrier to noise ratio in the 
loop) in (5.5.1). 


T ave " [2/ “n ]eX P 


(5.5.1) 
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The CNR^ was ca l cu l ate ^ using' the equivalent noise bandwidth 
(B^) for the linear baseband model of the PLL, This is derived by 
Viterbi [ 3 ] . 


B t - + a, /2 £]Hz . (5.5.2) 

J-i u u 

In order to achieve an effective approximation of white noise 

at the input to the PLL, o> of the loop was set at .1 and t was 

n 

set to .707. This gives an equivalent noise bandwidth of <333 
radians per second compared to a filter bandwidth of one radian per 
second . 

To compare the results of the model developed in this thesis to 
those of Sanneman and Rowbotham, the mean time to unlock is con- 
verted to an average cycle skipping rate. This rate is then used to 
predict the number of impulses in the time intervals over which the 
thesis model was run. This comparison is shown in figure 18. 

When comparing the results one must remember that the results 
from Sanneman and Rowbotham's work represent a mean time to unlock. 
Thus, if the loop skips two cycles in a relatively short time inter- 
val, this would correspond to unlocking only once. This bunching 
effect has been observed previously [ 1 2 ] and was observed with the 
model of this thesis . However, bunching does not begin to occur 
until several cycles are skipped. Further, the bunching effect 
increases as the CNR^ is decreased. Thus, for high values of 
CNR^ , one would expect the rates found by this model to agree 
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closely with that predicted using the mean time to unlock. At low 
values of CNR^ one expect the rate given by the thesis model 

to be higher than that predicted using the mean time to unlock. In 
addition, one would expect the difference to increase as CNK^ 
decreases. All these results are apparent in figure 18. 

A comment on the variation of the data is in order at this 

4 

point. Computer runs were made in segments of 10 seconds each. 

The number plotted represents the average number of impulses in four 
such runs . The number of PLL impulses observed on individual runs 
varied by as much as a factor of four. This is, in part, due to the 
fact that, for an 03 of one tenth, the decorrelation time of the 
loop is a significant fraction of the total simulation time. At 
high carrier to noise ratios this meant the number of impulses 
observed was not sufficient to obtain a reliable estimate of the 
impulse rate. At low carrier to noise ratios the bunching effect 
previously mentioned indicates all the impulses observed are not 
independent events. In this case, the run times were not long enough 
to give a number of independent events which would allow as accurate 
an estimate of the impulse rate as one might like. This accounts 
for the variance apparent in figure 18. Considering these factors, 
the agreement is very good. 

A second comparison which supports the model is shown in figure 
19 . This is a comparison of the average noise power out of the loop 



Impulses Out per 10 seconds 
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Carrier to Noise Ratio in the Loop (dB) 

PLL Impulses for u one tenth of IF Bandwidth 


figure 18 
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Output Noise Power vs. Carrier to Noise Ratio in the Loop 
(u)^ Is One Tenth of IF Bandwidth) 


figure 19 
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versus CNR^ with the same results obtained experimentally.* The 
experimental results are plotted with a constant offset. This is a 
result of differences in total noise powers. The significance of 
the curves is not in the exact values of output noise power , but in 
the fact that the shapes are almost the same. Further, the varia- 
tions appear at low values of CNR^ • I n this region accurate 
experimental measurements are very difficult to obtain. 

5 . 6 New Results 

This section presents results obtained with the computer model 

for (£> of the PLL equal to the IF bandwidth. Figure 20 is a plot 

of the cycle skipping rate versus carrier to noise ratio in the IF. 

The CNR was used as the base parameter because the loop noise 

bandwidth is a meaningful parameter only when the white noise 

approximation at the input is valid. The results here were taken 

4 

from two runs of 10 seconds for each point plotted. It was noted 

that the variation in the number of cycles skipped on each run was 

much smaller in this case than for the simulations run with co of 

n 

the PLL set to one tenth. The variation was generally on the order 
of one tenth the value shown. For a comparison, the PLL impulse 
rate predicted with Sanneman and Rowbotham’s mean time to unlock is 
also plotted. This assumes the carrier to noise ratio in the loop 

* The experimental results were obtained by A. Gilbert. They will be 
published in an upcoming report to the National Aeronautics and 
Space Administration under grant NGR-32-003-037 . A similar curve 
is also given in reference [ 1 3 ] . 
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IF Carrier to Noise Ratio (dB) 

PLL Impulse Rate for Equal IF Bandwidth 


figure 20 
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is equal to that in the IF. The rate of IF impulses observed, as 
defined by Rice, is also plotted on the same graph. 

Figure 21 shows the effect on the number of PLL impulses of 
varying £ with the IF bandwidth equal to a> n of the PLL. 

Figure 22 shows the variation in average noise power out of the PLL 
for these same parameters. These results are of a preliminary 
nature and are included to illustrate the types of information that 
may be obtained with this model . 



Average Noise Power Out PLL Impulses per 10 Seconds 
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Variation in PLL Impulse Rate 
with Damping Ratio (for u> 
Equal IF Bandwidth) n 

figure 21 



Damping Ratio 

Variation in Output Noise Power 
with Damping Ratio (for cj 
E qual IF Bandwidth) n 


figure 22 


CHAPTER VI 


CONCLUSION AND RECOMMENDATIONS 


6 . 1 Conclusion. 

The conclusion that may now be drawn is that a useful model of 
a second order PLL has been developed and verified. The model, as 
shown in this thesis, can provide data on performance of a second 
order PLL for a wide range of carrier to noise ratios and ratios of 
loop bandwidth to IF bandwidth. There are many areas in which this 
model is valid for which no analytical solutions predicting the 
performance of a second order PLL have been obtained. This is 
particularly important for the case where cd of the PLL is of the 
same order of magnitude as the IF bandwidth. 

6 . 2 Recommendations for Future Study 

Since the direct simulation approach used in this thesis has 
been shown to produce valid results, one area of future study would 
be a thorough study of the many numerical integration techniques to 
determine the method which will give the best trade off between 
machine computation time and accuracy. 

As mentioned in Chapter V, there was some output from the com- 
puter runs which was not analyzed for the purposes of this thesis . 
This includes the duration of IF impulses and the times of occurrence 
of IF impulses and PLL impulses. Analysis of this data may allow a 
determination of a general relationship between IF impulses and PLL 
impulses . 
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As mentioned in the Introduction, only minor modifications 

would be required to add modulation to the model. This would simply 

involve structuring the main portion of the simulation in the manner 

shown in Chapter III, The desired output could be calculated and 

subtracted from the actual output to leave a measure of noise and 

distortion. It should be noted that, the results presented in this 

thesis required approximately 15 hours of time on the IBM system 

4 

360 model 50, Each run of 10 seconds required from 4.5 to 5.5 
minutes of computer time. The additional computations required for 
modulation could easily double this figure. Thus^ detailed studies 
with modulation would require a faster computer. The problem is 
further complicated by the fact that the introduction of modulation 
increases the dimension of the parameter space of interest. 

The extension to higher order loops is also straightforward. 

One may change the lowpass filter in the PLL to model any linear 
filter of reasonable order. The only restrictions are the accuracy 
of the computer and the computer time required. 
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APPENDIX A 


COMPUTER PROGRAM FOR THE DETERMINISTIC MODEL 


DIMENSION A (1000) , B(1000) , Z(1000) C(1000) 

DO 100 JJ=1.,20 
READ(5,2) ZETA , S S , WN 
2 FORMAT (3F10. 5) 

WRITE (6, 40) ZETA , SS , WN 

40 FORMAT ( 1 1ZETA= f , FI 0 « 7 , 5X, ' SS= f , F10 .7, 5X, , WN= , ,F10.7) 
Cl =2 . 0*ZETA*WN 
C2=WN*WN 
T=0.0 
TH0=0 , 0 
THOl=0 .0 
THOO=0 o 0 
THOS0=0 „ 0 
THOS20=0.0 
DO 20 1=1,1000 
DO 10 J=l,5 
T=T+.05 
THI=SS*T 
THE=THI-TH01 
TH1=SIN(THE) 

THOSll=Cl*THl 

THOS 2 1-THOS 2 0+C 2* (THl+THO) *.025 
THOSKTHOS 1 1+THOS 2 1 
THO l.“THOO+(THOS 1+THOS 0) * . 025 
THO0=THOl 
THOS ORTHOS 1 
THOS 2 ORTHOS 21 
10 TH0=TH1 
Z ( I ) =T 
C (I ) -THOS 1 
A(I)=TH01 
20 B(I)<<CHE 

WRITE(6 , 50) (A(I),B(I),Z(I),C(I), 1=1, 1000) 

50 FORMAT (1H ,4E20.10) 

100 CONTINUE 
CALL EXIT 
END 



Block Diagram of Deterministic Model Annotated 
with Computer Program Variables 

figure 23 







APPENDIX B 


THE RANDOM NUMBER GENERATOR 

The program generating the gaussian numbers was obtained from 
the scientific subroutine package written for the computer used in 
this study „ The subroutine which produced the. normally distributed 
random numbers is called GAUSS. It utilizes another subroutine 
which generates uniformly distributed random numbers. This sub- 
routine is called RANDU. Listings of GAUSS and RANDU, along with 
short descriptions, can be found in the manual for this package [14]. 

The general approach of GAUSS is to combine 12 uniformly dis- 
tributed random numbers from RANDU to produce each normally distrib- 
uted random number , This is significant because the uniform random 
number generator is a pseudorandom sequence generator. That is, it 
generates a finite sequence of random numbers and then repeats the 
sequence. If a subsequence of length less than the repetition length 
of the main sequence is subjected to statistical tests it will have 

all the statistical properties of a sequence of random numbers. The 

29 

main sequence length with the subroutine used is 2 

Approximately 2 x 10~* gaussian random numbers were used for each 

£ 

noise sequence generated for use in the model. This means 2.4 x 10 

29 

numbers were used from RANDU. However, 2 is approximately 

g 

5.36 x 10 . Thus, more than 200 such sequences could be generated 
before the basic random number generator repeats . 



APPENDIX C 


THE NOISE GENERATOR 


DIMENSION A(1000) 

READ(5,1)IX 

1 FORMAT (110) 

B0=0. 1789234 
B1.=0. 7525181 
B2.=l .3095747 
B3=l .9373675 
B4=l . 1724909 
Y5I=0.0 

Y4 1=0.0 
Y3I=0 .0 
Y 21=0.0 
Y1I=0 . 0 
Y0l=0 .0 

DO 100 111=1,210 

DO 99 J=l,1000 

CALI GAUSS (IX, 1,0,0.0,X) 

Y50.=X-Y4I*B4 -Y3I*B3 -Y2l*B2 -Y1I*B1 -Y0I*B0 
Y40=Y4I+(Y504Y5I) * . 025 
Y30=Y3I+(Y40+Y4I) * . 025 
Y20=Y2I+(Y30+Y3I)*.025 
Y10=Y1I+(Y20+Y2I) * . 025 
Y00=Y0I+(Y10+Y1I) * . 025 
Y5I=Y50 
Y4I=Y40 
Y3I=Y30 
Y2I=Y20 
Y1I=Y10 
YOI=YOO 
99 A(J)=YOO 
WRITE (10) A 
IF (III .LT .10) GO TO 100 
100 CONTINUE 

WRITE (6, 2) IX 

2 FORMAT ( 1H ,120) 

STOP 

END 
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APPENDIX D 

RADIUS OF GYRATION OF NOISE FROM A FIFTH ORDER CHEBYSHEV FILTER 

The radius of gyration,, y } of a narrow-band noise process is 
given by: 



where G(o>) is the lowpass equivalent power spectrum of the narrow- 
band noise [ 15 ] . For convenience, the integral in the numerator will 
be called 1^ and the integral in the denominator will be called 
I 2 . The spectrum of the noise out of the Chebyshev filter used in 
this thesis is: 


G(cq) = tH 0 /2] ! H ( j o>) | 2 


(D.2) 


where H(jco) 
Now I 1 


is the transfer function of the Chebyshev filter, 
will be calculated: 


I 1 = N 0^ 2 / t jctHC jco) ] [- ja:H(-ja>) ]dco . (D.3) 


Substituting s for jao , (D.3) can be rewritten as: 


1 ± = [n q / 2] -2it 


CD 


2rtj / [sH(s) ] [-sH(-s)]ds 


'-CD 


(D.A) 
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The integral in brackets has been tabulated in a general form [16]. 
Using the values for a fifth order Chebyshev filter one can obtain, 

I t m (N 0 /2)-4Jt 2 (.0475) . (D.5) 

The integral is simply the equivalent noise bandwidth for a 

fifth order Chebyshev filter. This is also tabulated [ll]. From 

the tabulated value for = 1 . 0 ; 

c 

I 2 « 2(N 0 /2)(1.02) . (D. 6 ) 

Combining (D.5) and (D. 6 ) one obtains: 

7 = [2* (.0475/1.02)]^ . (D.7) 

This gives a 7 of .915, Rice’s work was done in terms of frequency 
in cycles as opposed to the radian frequency shown here. Thus, the 
r for Rice's formula, shown in (4.4.13), is 7 divided by 2 Jt . 
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APPENDIX E 


MAIN SIMULATION PROGRAM 


DIMENSION SS(1)_,XS (1000) SA( 1) _,XC (1000) 

SS(1)=0.0 
SA(1)=0 .0 
C 

C READ IN NUMBER OF SETS OF DATA 
G 

READ(5, 1)NDATA 

1 FORMAT (15) 

DO 100 NNN=1,NDATA 

READ INPUT DATA FOR RUN 

READ(5, 2)SN , WN , 2ETA 

2 FORMAT (3F1 0.5) 

WRITE OUT DATA FOR RUN 

WRITE ( 6 , 3) SN , WN, ZETA 

3 FORMAT ( 1 1 SN='F10.5,5X , WN= , F10.5,5X , ZETA='F10.5) 

COMPUTE GAIN TO GIVE DESIRED CARRIER TO NOISE RATIO 


ANP-0 . 0 
PI-3. 1415927 
PI02=PI/2.0 

GAIN=0 . 1789234*SQRT(20 .0*PI/(SN*2 ,0)) 

INITIALIZE VARIABLES 

NRI=0 

N«1 

NPLLP-0 
NPLLN=0 
NIMP=0 
NPLUS-0 
NMINUS=0 
MUE=0 
PI2=2 .0*PI 
AVINT-0 „0 
M=0 
T=0 „ 0 
THO 1=0.0 
THO0=0 . 0 
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THOS 0=0.0 
THOS .20=0.0 
TH0=0 . 0 

Cl.=2 , 0*ZETA*WN 
C2=WN*WN 


START MAIN LOOP 

DO 98 J=1 , 200 
READ(10)XS 
READ(11)XC 
DO 50 L=l, 1000 
XS(L)=XS(L)*GAIN 
50 XC (L ) =XC (L ) *GA IN 
DO 98 K=l,25 
DO 99 MM=1,40 
I.-MM+(K- 1 ) *40 
T~T+ . 05 
TH=-SIN (TH01) 

THl=TH-XS (I )*TH+XC ( I )*COS (THOl ) 

THOS 1 1=IH1*C1 

THOS 21 =THOS 2 0+ ( TH0+TH1 ) * . 025*C2 
THOS l=THOS 1 1+THOS 2 1 
THO 1 -THO 0+ ( T HOS 1+THOS 0)* . 025 
THO0=THOl 
THOS 0=THOS 1 
THOS20=THOS21 
TH0=TH1 
C 

C MAIN PORTION OF SIMULATION ENDS 
C BEGIN TESTS FOR CYCLE SKIPPING 
C 

ANP=ANP+THOS l*THOS 1*0 . 05 
AVINT=AVINT+TH01 
81 IF (XC(I) .LT. -1.000) GO TO 83 
IF (M 0 NE „ l)GO TO 99 
WRITER, 6) M,N,TST ; T 

6 FORMAT (3H M=, 15, 5X, 2HN=, 15, 5X,4HTST=,E15 . 7,5X,2HT=,E15 . 7) 
IF(N.EQ.l) GO TO 87 
NIMP=NIMP+1 

IF(XS(I).GT.0.0) GO TO 86 
WRITE (6, 200) 

200 FORMAT ( ' POSITIVE IF IMPULSE') 

NPLUS=NPLUS+1 
GO TO 87 
86 WRITE(6,300) 

300 FORMAT (' NEGATIVE IF IMPULSE') 

NMINUS=NMINUS+1 
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87 N=1 
H“0 

TE-T-TST 
WRITE ( 6/7) TE 

7 FORMAT (' IMPULSE DURATION® ^EIS .7) 

GO TO 99 

83 XF(XS(I)*XS(I-1) .GE.0.0) GO TO 84 
N~N:*(-1) 

NRT.«NRI+1. 

WRITE (6, 85 )T 

85 FORMAT (' OZERO CROSSINGS ' E15 »7) 

84 IF(M„NE a l) TST^T 
M-l 

99 CONTINUE 

AVINT®AVINT/40 . 

IF ( ('MUE+1 )*PI2 „GT .AVINT ) GO TO 90 
MITE -MUE+1. 

NPLLP =NPLLP+1 
WRITE (6,91)1, AVINT 

91 FORMAT ( 'OPOSITIVE PLL IMPULSE AT ' , 2E15 „7/) 

90 IF((MUE-l)*PI2.LT.AVINT)GO TO 98 

MUE=MUE-1 
NPLLN=NPLLN+1 
WRITE (6,92)T, AVINT 

92 FORMAT (' ONEGATIVE PLL IMPULSE AT',2E15.7/) 

98 AVINT-0 „0 

WRITE (6,4000) NPLUS , NMINUS , N IMP , NPLLP , NPLLN 

4000 FORMAT ( r NPLUS^= 1 , 15 , 5X / NMINUS™ ' , 15 , 5X, 1 NIMP= ' , 15 , 5X , ' NPLLP = r , 15 , 5 
IX, 'NPLLN™ 1 ,15) 

ANP-ANP / 10000 . 0 
WRITE (6, 4001) ANP 

4001 FORMAT (' ANP=',E15.7) 

WRITE (6, 4010) NR I 

4010 FORMAT (’ NUMBER OF RICE IMPULSES = ',I 10) 

WRITE (6, 9) NNN 

9 FORMAT (4H RUN, 12, 2X, 8HC0MPLETE) 

REWIND 10 
REWIND 11 
100 CONTINUE 
CALL EXIT 
END 



